使用Python Geotiff转换数字高程数据的基本地图信息 |
您所在的位置：网站首页 › python 安装elementtree › 使用Python Geotiff转换数字高程数据的基本地图信息 |
|
2020.06.10更新
由于此安排已从北向南颠倒，因此添加了紧急更正。 介绍 标题中的数据在此站点上分发。 提供日本非常详细，准确的海拔数据(DEM)。 内务和通信部统计局为每个二级分区网格下载了一个zip格式文件。它包含xml格式的高程数据。稍后将描述区域网格。 代码在本文结尾。一个例子如下。 (由于文件名也用作元信息，因此请将其应用到下载的文件中) 命令 1python zdem2tif.py FG-GML-1234-56-DEM10B.zip FG-GML-6543-21-DEM5A.zip ...每个zip文件(=次要分区区域网格)将输出一个geotiff图像。 支持5m和10m分辨率。 我们已经确认了python3系列的操作。必须安装gdal软件包。 这次，我对这个站点上的代码进行了很多引用。但是，如果保持原样，则在数据丢失的情况下会发生错误，因此该代码已得到修复。 似乎还有其他各种工具，但是创建此代码是因为它可以在命令上执行，可以与外国人共享，因为它是用英语编写的，因此无需解压缩。 另外，正如我在发布后注意到的那样，qiita已经有一篇有关基本地图信息和python的文章。文件格式的说明在此处更详细。 关于区域网格 本文档中给出了详细定义。 根据空间分辨率在1到3阶分区中定义区域网格。 主网格被分配了4位数字的分区代码，辅助网格被分配了另外2位数字，第三网格被分配了另外2位数字。 可以使用车厢代码的编号来知道车厢四个角的纬度和经度。 源代码说明 将来，分布式数据的规范可能会发生变化，并且有可能在作者未确认的领域中不受支持，因此我将简要解释代码。 (convert)的总流量为 从zip文件名获取区域网格信息和DEM分辨率 定义一个输出文件数组(到目前为止DEMOutArray，mesh_info_from_zipname) 关于zip文件(for filename in filelist)中xml文件的循环 从xml文件(XMLContents，read_xml)获取所需的信息 为每个xml文件(DEMInArray)生成DEM数组 覆盖输出数组(update_raster) 以Tiff格式保存(save_raster) 从zip文件名获取区域网格信息和DEM分辨率 读取4位主分区代码和2位辅分区代码以计算四个角的纬度和经度。 此外，通过将文件名写为DEM5X还是DEM10X来解释分辨率。 对于每种分辨率，一个xml文件的数组大小都是固定的(例如，在10mDEM的情况下，假定zip中只有一个xml)。 定义输出文件数组预先创建一个次级部分的大小。 在5mDEM的情况下，假设一个zip(第二分区)中最多包含10 x 10 xml(第三分区)。但是，似乎有许多区域尚未被观察到，并且在某些情况下仅存在约10个xml。 -9999的异常值用于对应于不存在xml的区域的位置。 关于zip文件中的xml文件Loop第一次，我了解了默认包含的模块 zipfile。 通过使用此功能，您无需解压缩每个文件就可以读取zip中的文件。 从xml文件获取所需的信息 对于 5mDEM，请读取文件名中描述的第三级分区代码的最后两位，并确定输出数组中的插入位置。如上所述，在10mDEM的情况下，首先假设没有第三级分区代码的描述，而只有一个xml。 暂时获取文本类型的必要信息。假定xml等中的标记规则如代码中所示。如果您使用不同的规则阅读xml，它应该立即显示为错误。 为每个xml文件生成DEM数组3??769113使用 xml文件中描述的信息生成数组。 文件规范具有DEM数据的定义。据此，每个像素的土地分割是{地表，表面，海平面，内陆水面，无数据等}，并且对于值未定义的地方输入-9999的数量。 。 在最初引用的代码中，根据是否为"地面"来获取DEM信息，但是在此代码中，根据xml中的DEM值是否为-9999对其进行划分("其他")。因为在指示的区域中有DEM值。说"其他"的原因未知)。 对于未定义DEM的地方，将应用离群值-1，并将离群值的使用与首先没有上述xml文件的情况分开。 覆盖输出阵列关于 数组的方向，假定纬度方向定义为北。从xml获得的数据中有gml:sequenceRule，并且将其写为'+x-y'的事实也支持它。如果这是'+x+y'等，则可能是朝南，但我尚未围绕该条件进行划分。 以Tiff格式保存无话可说。 输出tiff的图 转换后的文件，其顶部为5mDEM，底部为10mDEM，下载了相同的二级分割区域(几乎没有编辑，如轴值)。负值被掩盖。可以看出，每种分辨率下数据的丰度都因区域而异。 在下载之前很高兴知道，但是某个地方缺少任何信息吗？ .. ..
整个源代码 最后，整个代码。 zdem2tif.py 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197import os import sys import zipfile import numpy as np import xml.etree.ElementTree as et import osgeo.gdal import osgeo.osr class XMLContents(object): def __init__(self): self.ns = { 'ns': 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema', 'gml': 'http://www.opengis.net/gml/3.2', 'xsi': 'http://www.w3.org/2001/XMLSchema-instance', 'xlink': 'http://www.w3.org/1999/xlink' } def read_xml(self, zf, filename, dem_mode): if dem_mode == 'DEM10': self.posix = 0 self.posiy = 0 elif dem_mode == 'DEM5': _split = os.path.basename(filename).split('-') self.posix = int(_split[4][1]) self.posiy = int(_split[4][0]) with zf.open(filename, 'r') as file: self.text = file.read().decode('utf_8') self.root = et.fromstring(self.text) self.dem = self.root.find('ns:DEM', self.ns) self.coverage = self.dem.find('ns:coverage', self.ns) self.envelope = self.coverage.find( 'gml:boundedBy//gml:Envelope', self.ns) self.lower = self.envelope.find('gml:lowerCorner', self.ns).text self.upper = self.envelope.find('gml:upperCorner', self.ns).text self.grid = self.coverage.find( 'gml:gridDomain//gml:Grid//gml:limits//gml:GridEnvelope', self.ns) self.low = self.grid.find('gml:low', self.ns).text self.high = self.grid.find('gml:high', self.ns).text self.tuplelist = self.coverage.find( 'gml:rangeSet//gml:DataBlock//gml:tupleList', self.ns).text self.gridfunc = self.coverage.find( 'gml:coverageFunction//gml:GridFunction', self.ns) self.rule = self.gridfunc.find('gml:sequenceRule', self.ns) self.start = self.gridfunc.find('gml:startPoint', self.ns).text if self.rule.attrib['order'] != '+x-y': print('warning sequence order not +x-y') if self.rule.text != 'Linear': print('warning sequence rule not Linear') file.close() class DEMInArray(object): def __init__(self, contents): self._noValue = -1. self.llat, self.llon = np.array( contents.lower.split(' '), dtype=np.float64) self.ulat, self.ulon = np.array( contents.upper.split(' '), dtype=np.float64) self.lowx, self.lowy = np.array( contents.low.split(' '), dtype=np.int) self.highx, self.highy = np.array( contents.high.split(' '), dtype=np.int) self.sizex, self.sizey = self.get_size() self.x_init, self.y_init = np.array( contents.start.split(' '), dtype=np.int) self.datapoints = contents.tuplelist.splitlines() self.raster = self.get_raster() def get_size(self): sizex = self.highx - self.lowx + 1 sizey = self.highy - self.lowy + 1 return sizex, sizey def get_raster(self): _x, _y = self.x_init, self.y_init raster = np.zeros([self.sizey, self.sizex]) raster.fill(self._noValue) for dp in self.datapoints: if _y > self.highy: print('DEM datapoints are unexpectedly too long') break s = dp.split(',') if len(s) == 1: continue desc, value = s[0], np.float64(s[1]) # if desc == '地表面': raster[_y, _x] = value if value > -9998 else self._noValue _x += 1 if _x > self.highx: _x = 0 _y += 1 return raster class DEMOutArray(object): def __init__(self): self._fillValue = -9999. def initialize_raster(self): raster = np.zeros([self.sizey, self.sizex]) raster.fill(self._fillValue) return raster def get_trans(self): # +x-y return [self.llon, self.resx, 0, self.ulat, 0, -self.resy] def update_raster(self, tile, posiy, posix): xmin = posix * tile.shape[1] ymin = posiy * tile.shape[0] xmax = xmin + tile.shape[1] ymax = ymin + tile.shape[0] self.raster[ymin:ymax, xmin:xmax] = tile[::-1] def get_size(self): sizex = (self.highx - self.lowx + 1) * self.tilex sizey = (self.highy - self.lowy + 1) * self.tiley return sizex, sizey def mesh_info_from_zipname(self, zipname): ''' Ref: Table 4 of https://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf ''' _split = os.path.basename(zipname).split('-') _lat1 = int(_split[2][:2]) _lon1 = int(_split[2][2:]) _lat2 = int(_split[3][0]) _lon2 = int(_split[3][1]) self.llat = _lat1 * 0.66666666666 + _lat2 * 0.08333333333 self.llon = _lon1 + 100. + _lon2 * 0.125 self.ulat = self.llat + 0.08333333333 self.ulon = self.llon + 0.125 if _split[4][:5] == 'DEM10': self.mode = 'DEM10' self.tilex, self.tiley = 1, 1 self.lowx, self.lowy = 0, 0 self.highx, self.highy = 1124, 749 self.sizex, self.sizey = self.get_size() self.resx, self.resy = 1.1111111111e-04, 1.1111111111e-04 elif _split[4][:4] == 'DEM5': self.mode = 'DEM5' self.tilex, self.tiley = 10, 10 self.lowx, self.lowy = 0, 0 self.highx, self.highy = 224, 149 self.sizex, self.sizey = self.get_size() self.resx, self.resy = 5.5555555555e-05, 5.5555555555e-05 else: raise ValueError('Unexpected definition of DEM:', _split[4]) self.raster = self.initialize_raster() def save_raster(self, out_filename): trans = self.get_trans() srs = osgeo.osr.SpatialReference() srs.ImportFromEPSG(4612) driver = osgeo.gdal.GetDriverByName('GTiff') output = driver.Create(out_filename, self.sizex, self.sizey, 1, osgeo.gdal.GDT_Float32) # output.GetRasterBand(1).WriteArray(self.raster) |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |